Examen final Modules 4 et 5
DUBii 2021
Emmanuelle Fouilloux-Meugnier
08 avril, 2021
Consignes
Complétez ce document en remplissant les chunks vides pour écrire le code qui vous a permis de répondre à la question. Les réponses attendant un résultat chiffré ou une explication devront être insérés entre le balises html code. Par exemple pour répondre à la question suivante :
La bioinfo c'est : <code>MERVEILLEUX</code>.N’hésitez pas à commenter votre code, enrichier le rapport en y insérant des résultats ou des graphiques/images pour expliquer votre démarche. N’oubliez pas les bonnes pratiques pour une recherche reproductible ! Nous souhaitons à minima que l’analyse soit reproductible sur le cluster de l’IFB.
Introduction
Vous allez travailler sur des données de reséquençage d’un génome bactérien : Bacillus subtilis. Les données sont issues de cet article :
Analyses
Organisation de votre espace de travail
# Je cree differents repertoires afin de stocker les données au fur et à mesure des etapes
mkdir -p ~/Module4-5/EvaluationM4M5-main/DATA/FASTQ
mkdir -p ~/Module4-5/EvaluationM4M5-main/DATA/CLEANING
mkdir -p ~/Module4-5/EvaluationM4M5-main/DATA/MAPPING
mkdir -p ~/Module4-5/EvaluationM4M5-main/DATA/QC
# Je me deplace dans mon repertoire et valide la presence des dossiers
cd ~/Module4-5/EvaluationM4M5-main
tree ~/Module4-5/EvaluationM4M5-mainTéléchargement des données brutes
Récupérez les fichiers FASTQ issus du run SRR10390685 grâce à l’outil sra-tools [1]
#Tout d'abord reserver les ressources de calcul avec salloc
salloc --cpus-per-task=6 --mem=1G
#charger sra-tools et identifier la version
module avail sra
##sra-tools/2.10.0 sra-tools/2.10.3
module load sra-tools/2.10.3
#avec fasterq-dump recuperer les données puis les compresser
srun --cpus-per-task=6 fasterq-dump --split-files -p SRR10390685 --outdir ./DATA/FASTQ
srun gzip *.fastqVerifier les fichiers telechargés en lisant les 6 premieres lignes.
cd ~/DATA/FASTQ
zcat SRR10390685_1.fastq.gz | head -n 6
zcat SRR10390685_2.fastq.gz | head -n 6Combien de reads sont présents dans les fichiers R1 et R2 ?
module load bc/1.07.1
echo $(zcat SRR10390685_1.fastq.gz | wc -l) /4| bc
echo $(zcat SRR10390685_2.fastq.gz | wc -l) /4| bcLes fichiers FASTQ contiennent 7066055 reads.
Téléchargez le génome de référence de la souche ASM904v1 de Bacillus subtilis disponible à cette adresse
wget ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/009/045/GCF_000009045.1_ASM904v1/GCF_000009045.1_ASM904v1_genomic.fna.gzQuelle est la taille de ce génome ?
#Utiliser l'outil seqkit
module load seqkit/0.14.0
seqkit stats GCF_000009045.1_ASM904v1_genomic.fna.gzLa taille de ce génome est de 4215606 paires de bases.
Téléchargez l’annotation de la souche ASM904v1 de Bacillus subtilis disponible à cette adresse
wget ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/009/045/GCF_000009045.1_ASM904v1/GCF_000009045.1_ASM904v1_genomic.gff.gzCombien de gènes sont connus pour ce génome ?
zgrep ID=gene GCF_000009045.1_ASM904v1_genomic.gff.gz | wc4536 gènes sont recensés dans le fichier d’annotation.
Contrôle qualité
J’aimerai avoir quelques informations sur ces 2 fichiers fastq
module load seqkit
srun seqkit stats --threads 1 *.fastq.gz| file | format | type | num_seqs | sum_len | min_len | avg_len | max_len |
|---|---|---|---|---|---|---|---|
| SRR10390685_1.fastq.gz | FASTQ | DNA | 7,066,055 | 1,056,334,498 | 35 | 149.5 | 151 |
| SRR10390685_2.fastq.gz | FASTQ | DNA | 7,066,055 | 1,062,807,718 | 130 | 150.4 | 151 |
On voit que dans les 2 fichiers on a le meme nombre de sequences. En revanche la taille minimale dans le #1 est de 35 alors que dans le #2 elle est de 130
Lancez l’outil fastqc [2] dédié à l’analyse de la qualité des bases issues d’un séquençage haut-débit
cd ~/Module4-5/EvaluationM4M5-main/DATA
salloc --cpus-per-task=8 --mem=1G
module load fastqc/0.11.9
srun --cpus-per-task 8 fastqc FASTQ/SRR10390685_1.fastq.gz -o QC/ -t 8
srun --cpus-per-task 8 fastqc FASTQ/SRR10390685_2.fastq.gz -o QC/ -t 8La qualité des bases vous paraît-elle satisfaisante ? Pourquoi ?
- Oui
- Non
car le Phred score median est toujours superieur à 30 comme le montre le graphique Per base sequence quality dans les FastQC Report
Lien vers le rapports MULTIQC
Est-ce que les reads déposés ont subi une étape de nettoyage avant d’être déposés ? Pourquoi ?
- Oui
- Non
car les reads entre 1_fastq et le 2_fastq n’ont pas la meme longueur minimale (35 versus 130). Il n’y a plus d’adaptateur sur les reads de 1_fastq, mais il reste des Illumina Universal Adapter sur les reads de 2_fastq
Quelle est la profondeur de séquençage (calculée par rapport à la taille du génome de référence) ?
#La taille du génome de reference est de 4215606 paires de bases
#Le nombre de reads est de 7066055 *2
#P= (7066055*2)*150/4215606
#502XLa profondeur de séquençage est de : 500 X.
Nettoyage des reads
Vous voulez maintenant nettoyer un peu vos lectures. Choisissez les paramètres de fastp [3] qui vous semblent adéquats et justifiez-les.
module load fastp
#noter la version de fastp: fastp --version #fastp 0.20.0
srun --cpus-per-task 8 fastp --in1 FASTQ/SRR10390685_1.fastq.gz --in2 FASTQ/SRR10390685_2.fastq.gz --out1 CLEANING/SRR10390685_1.cleaned_filtered.fastq.gz --out2 CLEANING/SRR10390685_2.cleaned_filtered.fastq.gz --html CLEANING/fastp.html --thread 8 --n_base_limit 5 --cut_mean_quality 30 --cut_window_size 8 --length_required 100 --cut_tail --overrepresentation_analysis --json CLEANING/fastp.jsonLes paramètres suivants ont été choisis :
| Parametre | Valeur | Explication |
|---|---|---|
| n_base_limit | 5 | pour eliminer les sequences polyN surrepresentees dans le 1_fastq |
| cut_mean_quality | 30 | la mediane des Phred score est superieure à 30. Les reads sont de bonne qualité |
| cut_window_size | 8 | fenetre glissante de 8 bases pour le calcul moyen de la qualité |
| length_required | 100 | la longueur minimale des reads doit etre de 100 |
| cut tail | parametre pour couper ici la queue | |
| overrepresentation_analysis | pour eliminer les sequences surrepresentees |
Ces paramètres ont permis de conserver 6777048 reads pairés, soit une perte de 4.09% des reads bruts.
Alignement des reads sur le génome de référence
Maintenant, vous allez aligner ces reads nettoyés sur le génome de référence à l’aide de bwa [4] et samtools [5].
cd ./MAPPING/
module load bwa/0.7.17
#il faut tout d'abord indexer le genome de reference
srun bwa index GCF_000009045.1_ASM904v1_genomic.fna.gz #GCF_000009045.1_ASM904v1_genomic.fna.fai
#reserver un job allocation adapté
salloc --cpus-per-task=32 --mem=8G
#debuter le mapping à partir des reads filtres
module load bwa mem2/2.2.1
srun --cpus-per-task=32 bwa mem GCF_000009045.1_ASM904v1_genomic.fna.gz ../CLEANING/SRR10390685_1.cleaned_filtered.fastq.gz ../CLEANING/SRR10390685_2.cleaned_filtered.fastq.gz -t 32 > SRR10390685_on_GCF_000009045.1_ASM904v1.sam
#creer le fichier bam
module load samtools/1.10
srun --cpus-per-task=8 samtools view --threads 8 SRR10390685_on_GCF_000009045.1_ASM904v1.sam -b > SRR10390685_on_GCF_000009045.1_ASM904v1.bam
#trier le fichier bam
srun samtools sort SRR10390685_on_GCF_000009045.1_ASM904v1.bam -o SRR10390685_on_GCF_000009045.1_ASM904v1.sort.bam
#indexer le fichier bam trié
srun samtools index SRR10390685_on_GCF_000009045.1_ASM904v1.sort.bam
#Enfin je genere un multiqc global
cd ~/Module4-5/EvaluationM4M5-main/DATA/
module load multiqc
srun multiqc -d . -o .Combien de reads ne sont pas mappés ?
cd ./MAPPING
srun samtools flagstat
SRR10390685_on_GCF_000009045.1_ASM904v1.sort.bam > SRR10390685_on_GCF_000009045.1_ASM904v1.sort.bam.flagstat
srun samtools idxstats SRR10390685_on_GCF_000009045.1_ASM904v1.sort.bam > SRR10390685_on_GCF_000009045.1_ASM904v1.sort.bam.idxstats13571369 reads totaux - 12826829 reads mappés=744540
744540 reads ne sont pas mappés.
Croisement de données
Calculez le nombre de reads qui chevauchent avec au moins 50% de leur longueur le gène trmNF grâce à l’outil bedtools [6]:
cd ..
#Chercher le gene trmNF dans le fichier d'annotation du genome complet
zgrep trmNF GCF_000009045.1_ASM904v1_genomic.gff.gz > trmNF.gff3
#Croiser le fichier d'alignement de tous les reads avec la sequence de trmNF et ne recuperer que les reads avec 50% d'overlap
module load bedtools
srun bedtools intersect -a ./MAPPING/SRR10390685_on_GCF_000009045.1_ASM904v1.sort.bam -b trmNF.gff3 -f 0.50 > trmNF_reads.bam
# Tri et indexage des reads
srun samtools sort trmNF_reads.bam -o trmNF_reads.sort.bam
srun samtools index trmNF_reads.sort.bam
srun samtools flagstat trmNF_reads.bam > trmNF_reads.bam.flagstat
2801 reads chevauchent le gène d’intérêt.
Visualisation
Utilisez IGV [7] sous sa version en ligne pour visualiser les alignements sur le gène. Faites une capture d’écran du gène entier.
J’ai eu du mal à generer un fichier indexé .fai de mon genome de reference. Apres plusieurs recherches, orientées par les messages d’erreur je l’ai généré en faisant ceci:
cd ./MAPPING/
zcat GCF_000009045.1_ASM904v1_genomic.fna.gz | bgzip -c > Genome.fa.gz
samtools faidx Genome.fa.gzMa capture d’ecran d’IGV:
Attention vue partielle des reads mappés
ImageFinale
References
A work by Migale Bioinformatics Facility
https://migale.inrae.fr
Our two affiliations to cite us:
Université Paris-Saclay, INRAE, MaIAGE, 78350, Jouy-en-Josas, France
Université Paris-Saclay, INRAE, BioinfOmics, MIGALE bioinformatics facility, 78350, Jouy-en-Josas, France
